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I. INTRODUCTION 
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Department of Physics, Texas A&M University, College Station, TX 77843, USA 

Symplectic integrators evolve dynamical systems according to modified Hamiltonians whose error 
terms are also well-defined Hamiltonians. The error of the algorithm is the sum of each error 
Hamiltonian's perturbation on the exact solution. When symplectic integrators are applied to 
the Kepler problem, these error terms cause the orbit to precess. In this work, by developing a 
general method of computing the perihelion advance via the Laplace- Runge-Lenz vector even for non- 
separable Hamiltonians, I show that the precession error in symplectic integrators can be computed 
analytically. It is found that at each order, each paired error Hamiltonians cause the orbit to precess 
oppositely by exactly the same amount after each period. Hence, symplectic corrector, or process 
\ integrators, which have equal coefficients for these paired error terms, will have their precession 

■») ■ errors exactly cancel after each period. Thus the physics of symplectic integrators determines the 

' optimal algorithm for integrating long time periodic motions. 
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Numerical methods for solving physical problems are generally not expected to contain interesting physics. They 
are viewed as mere means, or recipes, of arriving at a needed numerical solution. This is because most numerical 
methods are based on matching Taylor series, whose error terms have little to do with physics. By contrast, symplectic 
£^ . integrators solve dynamical problems by approximating the original Hamiltonian by a modified Hamiltonian whose 
error terms are also well-defined Hamiltonians. In the past, these error terms are just formal entities destine to be 
eliminated by order-conditions, and are rarely studied in their own right. Here, we show that a comprehensive study 
of the error Hamiltonians in the Kepler problem gives insights into the working of symplectic integrators and makes 
CD manifest, ways of optimizing them. 

Symplectic Integrators (SI)i*^ despite their excellent conservation properties, are not immune from the fundamental 
phase error when solving the Kepler problem. While the energy error is periodic, the phase error can accumulate 
and grow linearly with tim o 4 i 5 i 6 . One manifestation of the phase error is the "perihelion advance" of the numerically 
Oh' computed elliptical orbit. This error is particularly pernicious when contemplating long time integration of periodic 
motions. No matter how small the initial time step, the orbital precession error can accumulate after each period and 
grow linearly without bound. 

In the Kepler problem, the energy error causes the length of Laplace- Runge-Lenz (LRL) vector to oscillate and the 
phase error causes the vector to rotated. While the energy error has been studied extensively, little is known about 
the phase error and its cause. This is reflected in the historical development of symplectic integrators; most early 
integrators are not well-tuned for the reduction of phase errors. For example, when solving the Kepler problem, the 
first fourth-order, Forest-Ruth 8 algorithm has a much larger precession error per period than the standard fourth- 
order Runge-Kutta algorithm^. Even the much improved McLachan integrator— has a larger precession error than 
that of Runge-KuttaiS. 

In this work, we present a detailed study of the precession error due to each error Hamiltonian (up to fourth 
order) on Kepler's orbit. Based on Sivardiere's method^ 1 of computing the rotation of the LRL vector, we develop 
a comprehensive treatment of perihelion advance due to any perturbing Hamiltonian, including non-separable ones. 
We show analytically that paired error terms of the form {T, Q} and {V, Q} rotate the LRL vector oppositely by 
exactly the same amount after each period. Here T and V are the kinetic and potential energy functions of the Kepler 
Hamiltonian, {A, B}'s are Poisson brackets, and Q's are higher order Poisson brackets of T and V. Algorithms 
with equal coefficients for these paired error terms would therefore have their precession errors precisely cancel after 
each period. This class of algorithm has been previously identified^ as symplectic correcto»i2ii4, or processi^iSiii 
algorithms. Symplectic corrector algorithms were originally derived for their computational efficiency; this work 
further identify them as a class of integrators with periodic precession errors. Thus the physical effects of these error 
Hamiltonians provide the needed insight for devising optimal integrators with periodic energy and phase conservation. 

For the Kepler problem, highly specialized algorithmsi&iS can be devised to exactly conserve energy and the rotation 
of the LRL vector. However, these algorithms do not limit the growth of the phase error in time. At a given time, the 
particle is at the wrong point of the trajectory, despite the fact it is constrained to move on the correct trajectory. 
Also, the phase errors in these algorithms are only second order in At and cannot be systematically improved to 
higher orders. This work solves the Kepler problem to fourth order in both the energy and the precession error and 
illustrates a general philosophy of allowing the physics of the problem to dictate the type of algorithm to be used for 
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its solution. 

In the next section, we will summarize needed results on the error structure of symplectic integrators. This is 
followed by Section III where we derive analytical expressions for the rotation angle of the LRL vector per period due 
to error Hamiltonians up to the fourth order. In this work, we systematize and generalize Sivadiere's methodii of 
computing orbital precession to include any angular-momentum conserving Hamiltonians, even non-separable ones. In 
Section IV, we numerical verify these theoretical predictions. In Section V, we derive second and fourth-order corrector 
algorithms with demonstrated periodic precession errors. Some conclusions and directions for future research are given 
in Section VI. 



II. ERROR HAMILTONIANS OF SYMPLECTIC INTEGRATORS 



Symplectic integrators for evolving the standard Hamiltonian 

Jf(q,p)=T(p) + V(q), with T(p) = , (2.1) 

can be derived^ by approximating the system's short-time evolution operator via a product of elemental evolution 
operators e eT and e £ v via 



c 



N 

ir +^«J]e^ f e^^, (2.2) 



where each Lie operators Q associated with variable Q acting on any other dynamical variable W is defined by the 
Poisson bracket 

QW = {W,Q}. (2.3) 

For a given set of factorization coefficients {tj, Vi}, the product on the RHS of [j2.2|l then produces a ordered sequence 
of displacements 

Pi(e) = c pi = pi - e — , 
oqi 

dT 

q i (e)=c £T q l =q l +e—, (2.4) 

OPi 

which defines the resulting algorithm. For a more detailed description, see Refi and RefiiS. For the study of time- 
reversible Hamiltonians, we will only consider time-reversible, symmetric factorization schemes such that either t\ = 
and Vi — VN—i+i, U+i = tzv— i+i> or vn = and Vi = vn-%, U = tx-i+i- (The use of asymmetric schemes to study 
time-reversible Hamiltonians may introduce unphysical and unnecessary distortionSi of the phase space at fintie At.) 
The product of operators in l|2.2|) can be combined by use of the Baker-Campbell-Hausdorff (BCH) formula to give 

N 



i=l 

where Ha is Hamiltonian operator of the algorithm. By the repeated use of (|2.3|l . one can deduce the Hamiltonian 
function Ha corresponding to the Lie operator Ha'- 

H a = e T T + e v V + e 2 (e T Tv{T 2 V} + e V Tv{VTV}) 

+ £ ( CTTTTV 

{TT 3 V} + {VT 3 V} 

+ &TTVTV 

{T (T V) 2 } + evTVTV {V(TV) 2 }) + ... , (2.6) 

where {TTV} = {T, {T, V}}, {T(TV) 2 } = {T, {T, {V, {T, V}}}} etc., are condensed Poisson bracket notations. This 
is the Hamiltonian function conserved by the algorithm. The error coefficients ey, ey, evTVTV, etc., are algorithm 
specific, calculable from knowing the {ti,Vi} coefficients^. In particular, 

N N 

&t = /Ji, e v = }^Vj. (2.7) 

i=l i=l 
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Thus all algorithms must have ex — 1 = &v m order to reproduce the original Hamiltonian. This will always be 
assumed. The Poisson brackets reflect properties of the original Hamiltonianii£ 

{TTV}} = Pi V ijPj , 

{VTV}} = -ViVi, 

{T (T V) 2 } = -2 Pi {V ikj V k + V ik V kj ) Pj , 

{V(TV) 2 } = ZViVijVj, 

{TT 3 V} = PiPjPkPiVijki , 

{VT 3 V} = -3p iPj V ijk V k . (2.8) 

To emphasize that these error terms are Hamiltonians, we will also denote Httv — {T, {T, V}}, Httttv — {TT 3 V}, 
etc.. 

For a central potential 



V(q) = V(r), (2.9) 



one can easily verify that 

Vi = v'ii, 



Vij = USij + (V" -U)hrj, (2.10) 

V ijk = U'iSijh + S jk ri + S ki rj) + (V'" - 3U')TiijT k , (2.11) 

V m = r-WfaSu + 5 jk 6 u + Skidji) + {V"" - W" + ir^U'^r^A 

+ (U" - r^U^^ijYkYi + Sjkhri + 4;r;f; + SuTjT k + 5jir k Ti + 4z? 4 fj) (2.12) 



where we have defined 



V'(r) 

U{r) = — y-L. (2.13) 



The forms (|2.10(1 - (|2.12(1 are arranged such that the derivatives are manifestly correct in one dimension. For the Kepler 
problem, where 



V(r) = --, (2.14) 
r 



the error Hamiltonians up to the fourth order are: 



Httv = r 3 (5 lj - 3h*j)PiPj , (2-15) 

Hvtv = -r-\ (2.16) 

Httvtv = 4r- 6 (S i:j - 6r<r j)p iPj , (2.17) 

H VTVTV = -4r- 7 , (2.18) 

35 

Httttv = -§r~ 5 (8 tj 8 k i - 10 5 ij v k v l + —r^jr^i )p iPjP kPi, (2.19) 

Hvtttv = 9r- 6 (% -BtifypiPj (2.20) 

Note that Httv, Httvtv, Hvtttv are all quadratic in p characterize by two numbers n and a, 

h(n,a) = r~ n (6 i3 - aTiVj)piPj- (2-21) 

The case of n = a will be shown to be specially simple. 

III. PERIHELION ADVANCES AS PERTURB ATI VE ERRORS 

The basic idea of Sivardiere's method^ of determining the precession of the Kepler orbit via the rotation of the 
LRL vector 



A = p x L — r , 



(3.1) 
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where f = r/r, is to extract the time derivative of A in the form of 

A = fi x A, (3.2) 
thereby identifying the precession angular frequency S~2, and obtain the precession angle over one period by integrating 



A8= [ Sl(t)dt, (3.3) 
Jo 



where P is the period. For our purpose, we will generalize Sivardiere's approach to treat arbitrary but angular- 
momentum conserving forces, including non-separable Hamiltonians. 
For any Hamiltonian which leaves L invariant, 



A = p x L + x (r x r). (3.4) 



For the Kepler Hamiltonian, 



Ho = ^P 2 - -, (3.5) 
2 r 

r = p, p = =^A = 0. (3.6) 

If l|3.5|l is perturbed by a central force of the form 

p = -V«(r) = /(r)r, (3.7) 

then one has 

A = -/(r)L x f . (3.8) 

Without lost of generality, we can always assume that the unperturbed A lies along the x-axis such that A = ei, 
whose length is the eccentricity e of the orbit. Thus we can cast l|3.8|l in the form 13. 2|) with 

n = -/(r)-cos(0)L, (3.9) 

and 

f2w 



(3.10) 



A6 = - (-f(r)r 2 )cos(9)d9, 
e Jo 

where we have used L = r 2 0. If /(r) can be expanded in inverse powers of r via 

-f(r)r 2 =J2Kr- n , (3.11 



where n = 0, 1, 2, etc., then by the use of 

1 1 



- = -(l + ecos0) with p = L 2 = a(l - e 2 ), (3.12) 
r p 



where a is the semi-major axis, one obtains the closed-form result 



P 

where we have defined 



Ae = E^«( e )> ( 3 - 13 ) 



C n {e) = - (1 + e cos 6) n cos 9d9. (3.14) 
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In table 1, we list the required integral C n (e) up to n = 8. Notice that for an inverse-square force, n = and AO = 0. 
By partial integration, it is easy to see that 

SJe)= [ (1 + e cos 6) n sin 2 (B)dS = —^—C n+1 (e). (3.15) 
Jo n + 1 

From this, one can also derive the following recursion relation for C n (e): 

(1 + -^)C n+1 = (2 + -)C n - (1 - e 2 )C n _!. (3.16) 
n+1 n 

For HvtVi corresponding to —f(r)r 2 — 4r~ 3 we have 

A9vtv = ^C 3 (e) = 1^(1 + ie 2 ). (3.17) 

For Hvtvtv 1 corresponding to —f(r)r 2 = 4 • 7r~ 6 , we have similarly, 

AOvtvtv = ^-C 6 (e) = 4 ' 7 6 6?r (l + ^e 2 + ^e 4 ). (3.18) 
P b P b 2 8 

The other perturbing Hamiltonians are not local potentials, but are non-separable Hamiltonians with angular- 
momentum conserving equations-of-motion, 

P = /(r,p)r + .g(r,p)(p • r)p, 

r = -S(r,p)(p-r)r + /i(r,p)p. (3.19) 



In this case, we have 



A = -/(r, p)L x f + g(r, p)(p • r)p x L — M£i£1l x f . (3.20) 



The last and the third term can be treated as discussed above. It is only necessary to expand —fr 2 and —h in inverse 
powers of r and invoke 13.13|) . The middle term requires further attention. We rewrite it as 

A = .g(r,p)(p-r)(A + r) (3.21) 

The first term above has the exact solution 



A(i)=exp|7 9 (r,p)(pT)d<|A(0), (3.22) 



which induces no rotation on A and can be ignored. For the second term, relative to L x f , f lags 90° behind, so that 
the corresponding $7 is given by 

n = S (r,p)(p.r)icoB(0- J)L, (3.23) 
e 2 



with 



1 f P 

A6 = - 5 (r,p)(p-r)sin(0)di. (3.24) 
e Jo 



In doing the time integration, one can use the unperturbed Kepler orbit, with p ■ r = rr and 



Hence, 



L = t S m(e)e. (3.25) 
r p 



i r 27T 

A6=- g(r,p)r 3 sm 2 (9)d9. (3.26) 
P Jo 



If g can be expanded in inverse power of r such that 



gr 3 = P™ r ~ 



then again we have the closed-form result 

P n a i \ P" C"+i( e ) 



9 ^ W ^ n+1 



For the quadratic Hamiltonian /i(n, a), we have equations-of- motion of the form H3.20(l with 

- fr 2 = -nr-" +1 p 2 + a(n + 2)r-"- 1 (p-r) 2 , 
gr 3 = 2ar- n+1 , 
-h = -2r- n . 

The precession angle from r 3 g and —h can be read off directly: 

i\Ug = la = I - — , 



A6 h = -2 



p" n p" 

C n (e) 
P n 



These two contributions exactly cancel if n = a. 

Since the time integration can be done along the unperturbed Kepler orbit, we can replace 



P 2 = ---, (p • r) 2 = P V - L 2 

r a 



and reduce — fr to only a function of : 



fr 2 = 2(a(n + 2) - n) r - n - -(a(n + 2) - n)r"™ +1 - a(n + 2)L 2 r~ n - 1 

a 



yielding 

1 

P 

By the use of recursion relation , this can be simplified to 

A0 f = — 

1 P" 

For a = n, we just have 



A9 f = [2(a{n + 2) - n)C n - (a{n + 2) - n)(l - e 2 )C„-i - a(n + 2)C n 



+ 



(1 - -(n + 2))C n + (a- n )^^C n+1 
n n + 1 



Combining results (|3.3U|) . I|3.31|l and (|3.36|) . for Httv (a = n = 3), we have 

AOttv = -4rC 3 (e), 

P 6 

which is the exact negative of AOvtv- For Httvtv (pi = n = 6), we have 

AOttvtv = — — -Ce(e)i 

which is the exact negative of AOvtvtv ■ 
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For HvtttVi n = 6 and a = 3, we have 

A9 VTTTV = 9(A6>/ + A9 g + A9 h ) 

= 9 



— (^3(76 — ^-Ct) H 1 — 2— | 



P 

9-4 
9 • 4 • 12tt 



7 

<"',;!' ) + ~C 7 (e) 



(-1 



(3.39) 



For Httttv i we have 



/r 2 = 9-5r~ 4 


P A 


- 14p 2 (p- 


r) 2 - 


f 21(p-r) 4 


gr 3 = 3 • 4 • 5 r 


-4 


7(p-r) 2 - 
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-h = 9-4r~ 5 


P 2 


-5(p-rf 







(3.40) 



By use of (|3.32l) . all can be expressed in terms of r, yielding correspondingly 

9-8-5 



A0 



/ 



P 

9-7-5 



A9„ 



A9 h = 



^i(^C 6 -4(l-e 2 )C 5 -5C 7 
P v 3 

^i(-8C 6 + 4(l~e 2 )C 5 + 5C 7 ) 



(4C 6 - 4(1 - e 2 )C 5 + (1 - e 2 ) 2 C 4 ) 
(3C 8 -8C 7 + 4(l-e 2 )C 6 ), 



(3.41) 



The repeated use of the recursion relation (|3.16H to eliminate all terms except Cq and C 7 simplifies the above to 



9 • 4 / 4 
A9 S = W (C 6 + -C 7 



A9 g = —^2C 6 ~-C 7 



(3.42) 



finally giving 



AOttttv = (A9f + A9 g + A9 h ) 



9 • 4 r 

"p6 



C 6 (e) + -C 7 (e) 



(3.43) 



which is the exact negative of A9vtttv ■ 



IV. NUMERICAL VERIFICATIONS 

By monitoring the rotation of the LRL vector of a given algorithm when solving the Kepler problem, one can 
directly check the analytical results of the last section. For this purpose, it is useful to employ algorithms with only 
a single error Hamiltonian. For example, the second order algorithm / 



7>(e) = exp(-eV) exp^sT) exp{\eV ) exp( 7 -ef ) exp(±eV) 
b 2 3 2 o 



(4.1) 



has modified Hamiltonian 23 



72 



Hvtv + 0{e 4 ). 



(4.2) 



8 



Algorithm II, obtained by interchanging f <-> V, has Hamiltonian 

Hi 1 = H + ^H TTV + 0(e 4 ). (4.3) 

By running both algorithms at smaller and smaller e, and dividing the rotation angle of the LRL vector after one 
period by e 2 /72 until convergence is seen, we can directly test the predicted result (|3.17ll . For starting values of 
r = (10, 0) and p = (0, 1/10), such that p — L 2 = 1 and e = 0.9, we have the theoretical result 

A9vtv = -A9ttv = 45.33318 . (4.4) 

Algorithm I at s = P/10000 with double precision gives 

A0j = -45.33157. (4.5) 

Algorithm II at the same step size produces 

Ma = -45.33316 . (4.6) 

Both are in excellent agreement with the theoretical value, including the sign. Each algorithm causes the LRL vector 
(and hence the orbit) to rotate differently in time, but at the end of the period, both algorithms have rotated the 
LRL vector by the same amount. This is shown in Fig. 1. 

To test Httttv an d HvtttVi we consider the following symmetric, fourth-order forward^ algorithm, 

T = ... exp(ew V r + s 3 u U) exp(etiT) exp(eviV + e 3 uiU) exp(et 2 f) exp(ev 2 V + e 3 u 2 U), (4.7) 

where we have only indicated operators from the center to the right and where 

ViV + e 2 Ui U (4.8) 

indicates that one should update the momentum by compute the force from the effective potentia l- 2 '^' 2 '-' 

Vl V + e 2 Ul {V, {T, V}} = Vi V - e 2 Ul (VV) 2 . (4.9) 

Here, U = {V, {T, V}} and has nothing to due with the function defined in Section II. For positive coefficients {ti} 
and {vi}, 

3 1 8 125 1 

tl= 10' t2= 5' V0 =27' Ul = 432' W2= 16' (4 ' 10) 

3121 1145 409 , t , 

Un = , Mi = , Ito = , (4.11) 

1710720 2737152' 1520640 v ; 

we have algorithm /// with Hamiltonian 

H* A " = H + ^^Hvtttv + 0(e 6 ). (4.12) 

This forward time-step algorithm with only a single fourth-order error term can be easily converted to a sixth-order 
forward algorithm— by solving Hvtttv directly as discuss below. For a different set of coefficients 

*i = ^. *2 = ~, «o = ^(4V3-3), ^ = ^(^3-3), « 2 = i(V3-l), (4.13) 

" =98k)( 943 - 461 ^' ^li^ 1 " 2 ^ " 2 = 84)( 617 - 344 ^ ^ 
we have algorithm IV with Hamiltonian 

H 1 / = H - j^(7 - 4V3)H TTTTV + 0(e 6 ). (4.15) 
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For the same initial condition as before, we have 

ABttttv = -A0VTTTV = 5933.72 . (4.16) 

For /// and IV, we increase e to avoid machine errors. Running both algorithms at e = T/5000 gives 

A9 IH = -5933.77 and A9 IV = -5933.68, (4.17) 

both are in excellent agreement with the predicted value (|4.16|) . The rotation of the LRL vector in time is given 
in Fig. 2. Despite the more complicated structure of the fourth-order Hamiltonians, the resulting rotations of the 
LRL vector are very similar to the second order case. The only discernable difference is that since the fourth-order 
Hamiltonians are more singular, the LRL vector rotates over a much narrower range near mid period. 

It has been shown in Ref^2 that for positive coefficients, it is not possible to have both cttttv and cvtttv vanish 
and hence not possible to isolate the error Hamiltonian Httvtv or Hvtvtv by itself. (Using negative coefficients 
would entail too many operators with only numerical, rather than analytical coefficients.) However, since the effects 
of Httttv an d Hvtttv have been verified, one can check the theoretical results for Httvtv and Hvtvtv in 
combination with Httvtv and Hvtvtv i n a general fourth-order algorithm. We will do this in the next section. 
For future reference, for the same initial condition, one has 

— AOttvtv — AOvtvtv — 1812.98 . (4.18) 

For the second and fourth-order algorithms considered in this section, the error coefficients evTV, &ttv and evTTTVy 
s-ttttv j are of opposite signs, resulting in algorithms which rotate the LRL vector in the same direction. This is not 
accidental, but a basic feature of forward symplectic algorithms to be discussed in the next section. 



V. SYMPLECTIC CORRECTOR ALGORITHMS 



A general second-order, time-reversible algorithm has modified Hamiltonian, 

Ha = H + e 2 (eTTvHTTV + &vtvH V tv) + 0(e 4 ). (5.1) 
For example, the velocity form of the Verlet algorithm 

T vv (e) = exp(i £ V0 exp(ef) exp(i £ T>) (5.2) 

has Cttv = 1/12 and evTV = 1/24. This allows us to immediately predict that when it is used to solve the Kepler 
problem, its precession angle per period, after being divided by e 2 , must be A#ttv/24 = —1.8888. This is illustrated 
in Fig. 3. In order to eliminate this second order precession error, one must devise algorithms with cttv = &vtv ■ 
This requirement is the same as for being a second order symplectic correctoni 3 ^, or processi^iSiii algorithm. 
More generally, a symplectic integrator T of order n is a corrector kernel algorithm if it is such that the similarity 
transformed algorithm STS^ 1 is of order n + 2, where S is the corrector or processor. This is possible only for T 
having equal error coefficients^ for each pair of error terms {T, Q} and {V, Q}. When corrector algorithms are applied 
to the Kepler problem, the precession error in each order would cancel after each period and both the energy and the 
precession error would be periodic in time. 

However, it is not easy to satisfy this second-order "correctablility" requirement of 

e-TTV = e V TV- (5.3) 

If either {ti} > or {vi} > 0, Chmi2 and Blanes-Casas^ have proved that it is not possible to have cttv = &vtv ■ 
Moreover, a recent theorem^ has precisely stipulated that for positive factorization coefficients, evTV and &ttv must 
be separated by a finite, calculable gap. If errv = 0, then bvtv < and if evTV = 0, then e-TTV > 0. However, it 
is easy to force evTV to equal evrv if Hvtv = {V, {T, V}} can be directly added to the potential as done in (|4.9|l . 
For example, the Takahashi-Imada (TI) integrator—, 

T T i = exp Q £ t) exp (eV + ^e 3 [V, [f , V]]j exp Qef) , (5.4) 

has Cttv = &VTV = —1/24 = —0.0416667. Its LRL rotation angle in solving the Kepler problem is shown in Fig. 3. 
The precession error, like that of the energy error, now returns to zero. If {U, Vi} are allowed to be negative, then the 
following corrector algorithm can also be used, 

Tnf = ■■■ exp(evoV) exp(etiT) exp(eviV) exp(e< 2 T), (5.5) 
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with 



1 



v = 



2 - 2V3' 



*2 = 2^0, 



Ul = <i 



(5.6) 



and e^TV = &vtv = —0.0470817. Its precession error is also shown in Fig. 3, denoted as the non-forward (NF) 
algorithm. Since its error coefficients are very close to that of TI, its behavior is also very similar. Note that this 
non-forward algorithm requires three force evaluations (the minimum necessary) which is not very efficient. For three 
force evaluations, one can have a fourth-order algorithm without any second-order errors. Omelyan^ has recently 
shown that the force gradient in general can be extrapolated with only one additional force evaluation. Thus the effort 
in computing the force gradient is the same as the force. We conclude from this discussion that the TI integrator is 
likely the optimal second-order algorithm for integrating Keplerian orbits with two force evaluations. 
For a fourth-order time-reversible algorithm, the modified Hamiltonian is 



H A = H + e 4 ( 

&TTTTV HttTTV + SVTTTV HyTTTV 
+ STTVTV HttVTV + SVTVTV HvTVTV ) + 0(£ 6 ). 



(5.7) 



By knowing the error coefficients Cttttv ,&vtttv t&ttvtv and evrvTVi the precession error of any fourth-order 
algorithm can be predicted. For example, the well known Forest-Ruth algorithm^ has the same form as (I5.5|) . but 
with coefficients 



1 



ti 



1 



Vl 



2-2V3' 



. 2 V3 



Vl, 



(5.8) 



error coefficients 



gttttv = -0.00041376, cvtttv = -0.00868165, 
errvTV = 0.00702660, e VT vTV = -0.02604494, 



(5.9) 



and precession error 



A6fR = {&TTTTV — evTTTv) A9tTTTV + {&VTVTV ~ GTTVTv)^^VTVTV 

= 49.0593- 59.9580, 
= -10.8987, 



(5.10) 



which is in good agreement with the observed error— of -10.8890 computed at e = P/10000. In contrast, the forward 
algorithm C-£ 



T c = ... exp(ev V + e 3 u U) exp(etiT) exp(eviV + e^uiU) exp(et 2 T), 



where 



has error coefficients 



1 

4' 



Vl 



1 1 1 

Un = , lii = 0. ti = — , <2 = — . 

192' 3' 6' 



(5.11) 



(5.12) 



7 



cttttv = 
ettvtv 
and a precession error of only 



-0.000135, e V TTTV 



1 



-0.000116, 



51840 ' 8640 

= -0.000304, evTVTV = — — = -0.000239, 



23040 



46080 



A8c — {CTTTTV — SvTTTv) AOtTTTV + ( e VTVTV — GTTVTv) ^8 yTVTV 

= -0.114462 + 0.118033, 
= 0.003570, 



(5.13) 



(5.14) 
(5.15) 



which is more than three order-of-magnitudes smaller. This theoretical value is again in excellent agreement with the 
algorithm's actual error of 0.003565 at e = P/10000. Algorithm C uses only one more force gradient (and therefore 
only one more force) than FR. We have previously demonstrated^ that algorithm C's precession error is smaller than 
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recent fourth-order symplectic integrator proposed by McLachan 9 ., Blanes and Moan (recommended in Refi^) and 
Omelyan, Mrylgod and Folk 29 ; 30 . 

For a fourth-order algorithm, the precession error will return exactly to zero only if algorithm is correctable with 

&TTTTV = BVTTTV (5.16) 

gttvtv = cvtvtv ■ (5-17) 

This partly explains why algorithm C is so much better than algorithm FR: its error coefficients are more nearly 
equal. However, its unusually small precession error is due also to the near cancellation of two distinct error types in 
(ETTII) . 

The equality (|5TT|) can be easily satisfied by redistributing the gradient term in C. For example, by changing only 

1 a 1 

un = (1 — a) and u\ = , (5.18) 

v 7 192 2 192 y ' 



with 



one changes only 



a=A (5.19) 



e TT VTV = e VTV TV = -^tt: = -0.000260 . (5.20) 
oo4U 



The precession error now goes up to 

/ = (&TTTTV ~ BVTTTv)A0TTTTV 

= -0.1144622. (5.21) 

While this is in excellent agreement with the observed value of —0.1144619 at e = P/10000, this is clearly not an 
improvement over algorithm C. Instead of forcing only errvTV — s-vtvtv, one can also choose 

9 4 A8 T TTTv(e) , . 
a = ; — 5.22 

10 15 A9 VT vTv(e) 

so that total precession error vanishes for given initial choice of the eccentricity e. For e = 0.9, we have 

a = 0.027225479 . (5.23) 

Numerically, the precession error of this tailored algorithm returns to A8 = —2.11 x 10~ 6 after one period. Since 
a = corresponds to algorithm C, this algorithm differs only slightly from C. However, the slight change is essential 
for forcing the precession error to zero. Its precession error is compared to that of C in Fig. 4. 

The above tailored algorithm is not a general algorithm because it requires a priori knowledge of the eccentricity 
of the orbit. For a general corrector algorithm, one must enforce (|5.16f) in addition to l|5.17[) . As in the second order 
case, the equality (|5.16|) cannot be satisfied for forward algorithms. One must therefore keep one of the two error 
Hamiltonians. We keep the simpler Hvtttv and generalize (|5.11|l to 

T c = ... exp(ev Q V + e 3 u U) exp^iT) exp(e«iF + e 3 MiC>) exp(et 2 f) cxp(e 5 u;iW / ), (5.24) 

where we have denoted simply, W — Hvtttv- The coefficient w\ is chosen to satisfy (|5.16|l . 

Since Hvtttv is non-separable, one must solve the general equation-of-motion 13.19f) implicitly. However, since 
this error term is of order e 4 and has a small coefficient u>i, any low order scheme with at most 1 iteration is sufficient. 
We used a second-order implicit midpoint schema 3 -. (A second order method is needed to preserve time-reversibility. 
However, at e — P/10000, the results are unchanged even with no iteration, or with the use of the naive Euler 
algorithm.) For algorithm C H5.12fl with (|5.19|) . we must have u>i = —1/103680. The resulting precession error indeed 
returns to zero, however its error near t = P/2 is « 0.1, which is unacceptably large. By use of the one-parameter 
family of algorithm 4ACB as described in Ref»iS, we have found the following, likely optimal, fourth-order symplectic 
corrector algorithm 4S, 
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, .29 a 29 455 1 

U0 = (1 - a W Ul= 2 4608' a= 1102' Wl = -86400- (5 ' 26) 

Its precession error is compared to that of C and C in Fig. 4. Algorithm 4S's precession error returns to 3.1 x 10~ 6 
after one period and is never more than 8.9 x 10~ 3 at any time. Its error coefficients are 

&TTTTV = evTTTV = = 0.0000347, 

Zaovv 

53 

&ttvtv = evrvTV = Aonnan = 0.0001211 . (5.27) 
437760 

The algorithm evolves in time pcrscrving the constancy of the modified Hamiltonian (|5.7() , 

H (t) + e A H 4 (t) = H o (0) + e i H 4 (Q) + 0(e 6 ), (5.28) 
where H 4 is the total fourth order error function. It can be extracted as 

Ht(t) - Hi(0) = lim l(ff o (0) - flo(*)). (5-29) 

The right-hand-side is plotted in Fig. 5. Algorithm C"s error is slight higher than than of C, while the maximum 
error of 4S is approximately three times smaller than that of C. For a more general class of fourth order forward or 
gradient algorithms other then 4ACB, see Refsi 29 i 30 ' 31 . 



VI. CONCLUSIONS AND DIRECTIONS FOR FUTURE RESEARCH 



When solving physical problems, symplectic integrators approximate the original Hamiltonian by a modified Hamil- 
tonian with a well-defined error structure. For time-reversible integrators, the error Hamiltonians come in pairs in the 
form of {T, Qi} and {V, Qi}. There is a clear separation between the mathematics of the algorithm, which fixes the 
error coefficients e-TQi and &vQii and the physics of the problem, which determines the error Hamiltonians {T, Qi} 
and {V,Qi}. In the past, when symplectic integrators are studied as numerical methods, only the error coefficients 
are analyzed so that they can be set to zero. Here, by a well-chosen example, we have shown that the physical effects 
of the error Hamiltonians determine how the error coefficients should be chosen. That is, the underlying physics of 
the problem determines the best algorithm for its own solution. 

For solving celestial mechanics problems dominated by Keplerian orbits, this work shows that the optimal integrators 
at each order are symplectic corrector algorithms. Unfortunately, for forward algorithms without any unphysical 
backward intermediate time steps, this cannot be implemented without including extra error Hamiltonians. In second 
order, it is easy to include Hyrv, which is just a local potential. In fourth order, Hvtttv is a non-separable 
Hamiltonian too cumbersome to be solved in general. One must find ways of including Hvtttv without solving it 
directly. 

The analytical results for the precession error are useful for verifying numerical calculations, however, it is a tedious 
way of proving the equality A9tq { = —A6yQ i . It should be possible to prove this equality without explicitly evaluating 
individual precession angles. 

We have shown in Ref^fi that the phase error in the harmonic oscillator vanishes when eTQ t = &VQi ■ It was simply 
not realized in that context that Htq { and HyQ i are also generating exactly opposite phase angles. From these 
two examples, maybe one can prove that for a general Hamiltonian with periodic orbits, only symplectic corrector 
algorithms can yield periodic errors for both the action and the angle variable. 

Finally, this work demonstrated that one must rethink the usual practice of minimizing the sum-of-square of the 
error coefficients as a mean of optimizing algorithms. The error Hamiltonians are not random; they come in pairs with 
opposite signs. The error coefficients should therefore be chosen to be pair- wise equal, i.e., one should seek optimal 
algorithms within the class of corrector algorithms. 
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FIG. 1: The rotation of the Laplace-Runge-Lenz vector due to second order error Hamiltonian —Hvtv and Httv in algorithms 
I and II. Each algorithm rotates the LRL vector differently in time, but by exactly the same amount after one period. Most 
of the rotation takes place near the mid period. The solid line gives the theoretical value of -45.33318. 



TABLE I: Explicit expressions for the function C„(e) 



n 


C„(e) 








1 


7T 


2 


2tt 


3 


37T(l + i e 2 ) 


4 


47r(l + |e 2 ) 


5 


5tt(1 + |e 2 + ie 4 ) 


6 


6tt(1 + |e 2 + |e 4 ) 


7 


7^(l+ 1 # e 2 + fe 4 + ^) 


8 


87 r(l + 2ie 2 + fe 4 + i|) 
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FIG. 2: The rotation of the Laplace-Runge-Lenz vector due to fourth-order error Hamiltonians Hvtttv and —Httttv ■ 
Because the fourth-order error terms are more singular, the rotation takes place over a narrower range near mid period. The 
solid line gives the theoretical value of -5933.72 . 
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FIG. 3: The rotation of the Laplace- Runge-Lenz vector for three second-order symplectic algorithms: velocity- Verlet (VV), 
Takahashi-Imada (TI) and the non- forward corrrector algorithm (NF). The solid line gives the theoretical rotation value of the 
VV algorithm: A0 T tv/24 = -1.8888 . 
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FIG. 4: The rotation of the Laplace-Runge-Lenz vector for three fourth-order integrators: algorithm C, algorithm C ' with 
added gradient term to force the rotation angle back to zero, and the true symplectic corrector algorithm 4S. As with most 
integrators, algorithm C's precession error does not return to zero. 




FIG. 5: The energy error functions of algorithms C, C and 4S. Algorithm 4S's maximium error is three times smaller than 
that of C. 



